Collapse of differentially rotating neutron stars and cosmic censorship 
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We present new results on the dynamics and gravitational-wave emission from the collapse of differentially 
rotating neutron stars. We have considered a number of polytropic stellar models having different values of 
the dimensionless angular momentum J jNP , where J and M are the asymptotic angular momentum and mass 
of the star, respectively. For neutron stars with J /M^ < 1, i.e., "sub-Kerr" models, we were able to find 
models that are dynamically unstable and that collapse promptly to a rotating black hole. Both the dynamics 
of the collapse and the consequent emission of gravitational waves resemble those seen for uniformly rotating 
stars, although with an overall decrease in the efficiency of gravitational-wave emission. For stellar models with 
J/M^ > 1, i.e. "supra-Kerr" models, on the other hand, we were not able to find models that are dynamically 
unstable and all of the computed supra-Kerr models were found to be far from the stability threshold. For 
these models a gravitational collapse is possible only after a very severe and artificial reduction of the pressure, 
which then leads to a torus developing nonaxisymmetric instabilities and eventually contracting to a stable 
axisymmetric stellar configuration. While this does not exclude the possibility that a naked singularity can be 
produced by the collapse of a differentially rotating star, it also suggests that cosmic censorship is not violated 
and that generic conditions for a supra-Kerr progenitor do not lead to a naked singularity. 

PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Lz, 97.60.Jd 



I. INTRODUCTION 

Differentially rotating neutron stars, either dynamically sta- 
ble or unstable, can be the result of several astrophysical 
scenarios such as stellar-core collapse or binary neutron star 
mergers (see, e.g. Refs. lfll- [Tlll for some recent works). Be- 
cause of their differential rotation, these stars can support 
masses higher than if they were uniformly rotating lfl2l . Fur- 
thermore, if J and M are, respectively, the angular momen- 
tum and the gravitational mass of the star, differentially ro- 
tating models can reach values of the dimensionless spin pa- 
rameter J/A/^ > 1; this is not possible for stars in uniform 
rotation, at least when described by a hadronic equation of 
state (EOS), in which case J/AP < 0.7 lO. Finally, when 
hypermassive (that is, having masses larger than the one of 
the associated mass-shedding configuration in uniform rota- 
tion), differentially rotating stars can be related to events such 
as gamma-ray bursts (GRBs). Hypermassive differentially ro- 
tating neutron stars can indeed be formed after the merger of 
binary neutron stars and their eventual collapse can produce 
a spinning black hole surrounded by a hot and massive disk, 
which could generate the relativistic jets that are observed in 
short GRBs 111 . 

The investigation of the collapse of uniformly rotating neu- 
tron stars to rotating black holes was first studied in detail in 
Ref. 13, where a specific set of dynamically unstable mod- 
els was constructed (named Dl to D4) for a polytropic index 
N = 1. The region of dynamical instability to axisymmetric 



perturbations was probed by following the stability criterion 
of Friedman, Ipser, and Sorkin [15], who suggested that the 
turning point along constant angular-momentum sequences 
also marks the onset of a secular instability. Models Dl to 
D4 in Ref. ifTill were then chosen to be near the turning-point 
line, but with somewhat larger central densities so as to en- 
sure dynamical (and not just secular) instability (see Table 1 
and Fig. 1 of Ref. H)'. 

In a series of works starting with Ref. iflTll . the study of the 
collapse of differentially rotating neutron stars was initiated, 
both in axisymmetry and in three spatial dimensions (3D). In 
Ref. ITtIi . in particular, three different models were consid- 
ered, two of which were "sub-Kerr", i.e., with J/AP < 1 
and one was instead "supra-Kerr", i.e., with J/AP > 1. The 
collapse of the latter was obtained by artificially depleting the 
pressure by 99%, but no investigation was made on the sta- 
bility of the progenitor models. The results presented in iHtIi 
were based on the use of a nonconservative numerical scheme 
and of a polytropic EOS preventing the formation of strong 
shocks, which cannot be handled well by nonconservative nu- 
merical methods. In a subsequent work [|l8||, the collapse of 



A recent investigation by Takami, Rezzolla and Yoshida jlfll aimed at de- 
teraiining the neutral point along these sequences (i.e., where the eigen- 
frequency of the fundamental mode goes to zero), shows however that the 
secular stability does not coincide with the turning point, but must be at 
smaller central densities. 
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similar models was studied in the presence of viscosity and the 
first results were published, but only in axisymmetry, for the 
collapse of differentially rot ating ne utron stars in the presence 
of a poloidal magnetic field lfl9l - l23ll . 

The main goal of this work is to reconsider the effect of dif- 
ferential rotation on the collapse of rotating neutron stars and 
to present a systematic investigation of the stability properties 
of differentially rotating models. Special attention has been 
paid to the following two questions. (1) Do stellar models ex- 
ist with J/Af ^ > 1 and that are dynamically unstable to the 
collapse to black hole? (2) If a stellar model with J/A/^ > 1 
is induced to collapse, does it lead to a naked singularity, thus 
violating cosmic censorship? 

Overall, our results can be summarized as follows: 

f i) differentially rotating, sub-Kerr neutron-star models 
can be found that are dynamically unstable to the grav- 
itational collapse to rotating black holes; 

( ii) the dynamics of sub-Kerr models is very similar to 
the one of uniformly rotating models: nonaxisymmetric 
instabilities do not have sufficient time to grow and the 
stars collapse promptly to black holes; 

( Hi) the efficiency in the emission of gravitational radi- 
ation is comparable to that of uniformly rotating stars, 
although generically smaller since the collapse is gen- 
erally slower; 

(iv) differentially rotating, supra- Kerr models that are 
dynamically unstable to the gravitational collapse to ro- 
tating black holes could not be found; rather, all of the 
supra-Kerr models studied were found to be dynami- 
cally stable; 

(v) a supra-Kerr model can be induced to collapse only 
through a very severe depletion of the pressure support 
which, however, does not lead to the prompt formation 
of a rotating black hole; 

f vi) because of the development of nonaxisymmetric in- 
stabilities, the gravitational-wave emission from supra- 
Kerr models could in principle be considerably larger 
than that of sub-Kerr models if ever induced to collapse. 

The paper is organized as follows. In Sec. we briefly 
summarize the equations solved and the numerical infrastruc- 
ture used. In Sec.|III]we present our study on the stability of 
differentially rotating neutron stars with different degrees of 
differential rotation and different polytropic EOSs. In Sec.lIVI 
we present the dynamics of the collapse of four models, three 
wifli J/M^ < 1 and one with J/AP > 1, and in Sec.|V]their 
gravitational-wave signal. Finally, in Sec.[VT]we summarize 
and conclude. 

Throughout this paper we use a spacelike signature of 
(— , +, +, +) and a system of units in which c ^ G = Mq = 
1. Greek indices are taken to run from to 3, Latin indices 
from 1 to 3 and we adopt the standard convention for the sum- 
mation over repeated indices. 



II. BASIC EQUATIONS AND NUMERICAL SETUP 

All the simulations presented here were done using the 
Whisky code which solves the general-relativistic hydro- 
dynamic equations on a three-dimensional numerical grid 
with Cartesian coordinates ll24ll . The code has been con- 
structed within the framework of the Cactus Computational 
Toolkit ll25ll . which provides high-level facilities such as par- 
allelization, input/output, portability on different platforms 
and several evolution schemes to solve general systems of par- 
tial differential equations. Clearly, special attention is dedi- 
cated to the solution of the Einstein equations, whose matter 
terms in nonvacuum spacetimes are handled by the Whisky 
code. 

In other words, while the Cactus code provides at each 
time step and on a spatial hypersurface the solution of the Ein- 
stein equations 

G^u ~ SttT^i, , (1) 

where G^i, is the Einstein tensor and is the stress-energy 
tensor, the Whisky code provides the time evolution of the 
hydrodynamic equations, expressed through the conservation 
equations for the stress-energy tensor T'^'^ and for the matter 
current density J'', 

y^Tf^" = , = . (2) 

This system of equations is then closed by an EOS which re- 
lates the pressure to the rest-mass density and to the specific 
internal energy. 

In what follows, and mostly for the sake of complete- 
ness, we give a brief overview of how both the right and the 
left-hand side of Eq. ([T]i are computed within the coupled 
Cactus /Whisky codes. The equations presented have al- 
ready been discussed in several different publications, e.g. in 
lITii 1261 127I1 and we refer the interested reader to these works 
for more details. We note that the Whisky code can also 
solve the equations of general-relativistic magnetohydrody- 
namics (GRMHD) within the ideal-MHD Umit (Sill Hi. 
Hereafter, however, in order to build the necessary under- 
standing of the dynamics of gravitational collapse in the pres- 
ence of differential rotation, we will consider unmagnetized 
fluids only, leaving the inclusion of magnetic fields to a future 
study. 

A. Evolution of the field equations 

We use the conformal and traceless decomposition of the 
Arnowitt-Deser-Misner (ADM) formulation of the Ein- 
stein equations as first presented in 3D in Ref. ll30ll . which 
is based on the ADM construction and has been further de- 
veloped in llsill . Details of our particular implementation of 
the conformal traceless reformulation of the ADM system as 
proposed by ll3^ - {32ll are extensively described in |j26l HH] and 
will not be repeated here. 

The code is designed to handle arbitrary shift and lapse con- 
ditions, which can be chosen as appropriate for a given space- 
time simulation. More information about the possible families 
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of spacetime slicings which have been tested and used with 
the present code can be found in ll26l[34ll . Here, we Umit our- 
selves to recalHng details about the specific foliations used in 
the present evolutions. In particular, we have used hyperbolic 
iC-driver slicing conditions of the form 

dta = -f{a)a^{K-Ko), (3) 

with f{a) > and Kq = K{t = 0), with K being the trace of 
the extrinsic curvature and a the lapse function. All the sim- 
ulations discussed in this paper were performed using condi- 
tion (O with / = 2 /a. For the spatial gauge we use one 
of the "Gamma-driver" shift conditions proposed in ll34ll (see 
also iHI). In particular, all the results reported here have been 
obtained using the hyperbolic Gamma-driver condition, 

dfi3' = Fdtr-ydti3\ (4) 

where /?' is the shift, are the "conformal connection func- 
tions" and F and rj are, in general, positive functions of space 
and time. For the hyperbolic Gamma-driver conditions it is 
crucial to add a dissipation term with coefficient rj to avoid 
strong oscillations in the shift. Experience has shown that 
by tuning the value of this coefficient it is possible to almost 
freeze the evolution of the system at late times. We typically 
choose F ^ 3/4: and r/ = 3 and do not vary them in time. 

The singularity-avoiding properties of the above gauge 
choices have proved equally good both when using excision, 
as it was done in Refs. Ill4ll and 12711 . and when not using 
excision jssll . In particular in this paper we employ the "no- 
excision" technique introduced in Ref. (251 and we add an ar- 
tificial dissipation of the Kreiss-Oliger type P6J on the right- 
hand sides of the evolution equations for the field variables (no 
dissipation is introduced for the hydrodynamical variables). 
As first pointed out in Ref. ll35ll . in fact, renouncing to excision 
and using instead suitable "singularity-avoiding" slicing con- 
ditions improves dramatically the long-term stability of the 
simulations, allowing for the calculation of the gravitational 
waveforms well beyond the quasi-normal-mode ringing. 

B. Evolution of the hydrodynamic equations 

An important feature of the Whisky code is the imple- 
mentation of a conservative formulation of the hydrodynamic 
equations 13711 . in which the set of Eqs. (|2]i is written in a 
hyperbolic, first-order, and flux-conservative form of the type 

9tq + a,f«(q) =s(q) , (5) 

where f (q) and s(q) are the flux vectors and source terms, 
respectively ifssll . Note that the right-hand side (the source 
terms) does not depend on derivatives of the stress-energy ten- 
sor 

An important feature of this formulation is that it allows for 
the extension to a general-relativistic context of the powerful 
numerical methods developed in classical hydrodynamics, in 
particular, hig h-resolution shock-capturing schemes based on 
exact lf39l - l4l[l or approximate Riemann solvers (see Ref. llssll 



for a detailed bibliography). Such schemes are essential for a 
correct representation of shocks, whose presence is expected 
in several astrophysical scenarios. 

For all the results presented here, we have solved the hydro- 
dynamic equations emp loying the Marquina flux formula and 
a third-order PPM ll42ll reconstruction. A third-order Runge- 
Kutta scheme was then used for the evolution. 

C. Mesh Refinement 

We solve both the spacetime and hydrodynamic equations 
on nonuniform grids using a "box-in-box" mesh refinement 
strategy implemented in Whi sky via the Carpet driver ll43ll . 
This introduces two important advantages: first, it reduces 
the influence of inaccurate boundary conditions at the outer 
boundaries, which can be moved far from the central source; 
second, it allows for the wave zone to be included in the com- 
putational domain and thus for the extraction of important 
information about the gravitational-wave emission produced 
during the collapse. In practice, we have adopted a Berger- 
Oliger prescription for the refinement of meshes on different 
levels Il44ll and used the numerical infrastructure described 
iniH. 

D. Gravitational-wave extraction 

While several different methods are possible for the extrac- 
tion of the gravitational-radiation content in numerical space- 
times, we have adopted a gauge-invariant approach in which 
the spacetime is matched with the nonspherical perturbations 
of a Schwarzschild black hole (see Refs. Ii45i , AQ for appli- 
cations to Cartesian coordinates). In practice, a set of "ob- 
servers" is placed on 2-spheres of fixed coordinate radius Tox, 
where the gauge-invariant odd-parity and even-parity 

^'^^^ metric perturbations (^Tj-Bst are extracted. From these 
quantities it is possible to compute the even- and odd-parity 
perturbations 

QL^^'^fl^ (6) 

where A = ^/2{£ + 2)!/(£ - 2)!. Using these quantities it is 
also possible to compute the gravitational-wave amplitudes in 
the two polarizations h+ and /ix as 

'^+-i^x -^Efe«-i fQUt')dA -X"^. (8) 

where ^y^™ are the s = — 2 spin- weighted spherical har- 
monics. We note that the approach discussed above cannot be 
employed when simulating the supra-Kerr model, for which a 
very high resolution is necessary and hence the use of a com- 
putational domain with rather close boundaries. As we will 
comment later on in Sec. IV Bl for that model the gravitational- 
wave emission will be computed through the quadrupole for- 
mula. 
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III. INITIAL STELLAR MODELS 

A. Equation of state and differential-rotation law 

We construct our initial stellar models using the numerical 
code rns [50] (see [51,^52[l for evaluations of its accuracy) as 
isentropic, differentially rotating relativistic polytropes, satis- 
fying the polytropic EOS 



p = Kp 



P + 



r - 1' 



(9) 



(10) 



where p, p, and e are pressure, rest-mass density, and energy 
density, respectively, while K = 100 and T are the polytropic 
constant and the polytropic exponent, respectively. In our 
discussion we will also use the polytropic index TV, defined 
through the relation T ee 1 + 1 /N. 

We further assume that the equilibrium models are station- 
ary and axisymmetric, so that the spacetime geometry is de- 
scribed by a metric of the form 



D^^^ (d0- wdt)2 +c2^(dr2 +r2d6|2), (11) 



where p and uj are functions of the quasi-isotropic co- 
ordinates r and 6 only. The degree of differential rotation, 
as well as its variation within the star, are essentially un- 
known and because of this we here employ the well-known 
"y'-constant" law of differential rotation Ii53i] 



{n — Ld)e 



2ip 



1 - {n-u})e 



2Tp ■ 



(12) 



where Vl = fl{r,9) is the angular velocity, ft^ is the angu- 
lar velocity at the center of the star, and ^ is a constant (with 
the dimensions of length) that represents the lengthscale over 
which the angular velocity changes (see ifsill for a new and 
more realistic law of differential rotation). In the remainder 
of this Section, we will measure the degree of differential ro- 
tation by the rescaled quantity A = A/r^, where is the 
equatorial coordinate radius of the star. For A — > oo uniform 
rotation is recovered while a low value of A indicates a high 
degree of differential rotation. 



B. Equilibrium models 

As mentioned in the introduction, when studying the gravi- 
tational collapse of a rotating neutron star to a Kerr black hole, 
an interesting question is about what happens to a configura- 
tion with an initial dimensionless spin parameter J/AP > 1 
{i.e., a supra-Kerr model). If the cosmic censorship conjecture 
is expected to hold, such models cannot collapse promptly to 
Kerr black holes, which are limited to J/A'P < 1. Rather, 
one expects that the transition to a black hole, when it occurs, 
takes place after the stellar model has shed some of its angu- 
lar momentum. This was indeed what was shown to happen 



in the previous study of Ref. lllTll . where a dynamically sta- 
ble supra-Kerr model was induced to collapse after a dramatic 
depletion of the pressure support. 

Here, we consider again this interesting question but, be- 
fore discussing in more detail the features of the gravita- 
tional collapse of either supra- or sub-Kerr models, it is 
worth discussing the properties of the equilibrium models 
that can be computed with the prescriptions of the EOS and 
of the differential-rotation law discussed in the previous sec- 
tion llll Al As will become apparent, the study of these equilib- 
rium models will be quite revealing for the stability properties 
and hence for what is realistic to expect from the gravitational 
collapse of differentially rotating neutron stars. To this end, 
we have constructed a large set of initial models for various 
values of the polytropic index N and degree of differential ro- 
tation A, reaching close to the mass-shedding limit and span- 
ning a wide range of central densities. 

Figure [U shows the value of J/AP as a function of cen- 
tral rest-mass density pc for the three different EOSs with 
N = 0.5, N 0.75, and N = 1.0. In these sequences 
the rotation law and the polar-to-equatorial coordinate axis 
ratio are fixed to A = 1.0 and rp/r^ = 0.35, respectively. 
The choice of A = 1.0 is a typical one, representing moder- 
ate differential rotation (the angular velocity at the axis and 
at the equator differ by a factor of ^ 3), while the axis ratio 
of 0.35 refers to very rapidly rotating models near the mass- 
shedding limit (when the limit exits). Along each sequence, 
we mark by a filled circle the model which roughly separates 
stable models (at lower central rest-mass densities) from un- 
stable models (at higher central rest-mass densities). As we 
do not know precisely which are the marginally stable models 
(no simple turning-point criterion exists in the case of differ- 
ential rotation), we use as a reference the stability limit of the 
nonrotating models and thus we mark with a circle the central 
rest-mass density of the nonrotating model having the maxi- 
mum mass for each EOS^. Stated differently, all models to the 
right of the circles are expected to be dynamically unstable or 
at least very close to the instability threshold. 

As becomes clear from this figure, all unstable models we 
were able to construct have J/AP < 1, i.e., are sub-Kerr. In 
contrast, in order to find supra-Kerr models, one must reach 
very low central rest-mass densities, where equilibrium mod- 
els are instead expected to be (very) stable against axisymmet- 
ric perturbations. Interestingly, for the particular sequences 
considered here, the value of J jAP in the unstable region 
{i.e., to the right of the filled circles) becomes nearly constant 
for all the considered values of the polytropic indices consid- 
ered; this represents an additional evidence that all unstable 
models are indeed sub-Kerr. 



^ We also note that the value of the central rest-mass density separating sta- 
ble from unstable models in the case of uniformly rotating neutron stars 
changes by less than K. 5% when moving from a nonrotating sequence to 
a maximally rotating one. As a result, using a nonrotating model to mark 
the stability is a very good approximation which does not affect the results 
presented here. 
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FIG. 1: J/M^ as a function of central rest-mass density pc for mod- 
els with polytropic indices = 0.5, A'^ = 0.75 and N = 1.0, when 
the rotation law and the polar-to-equatorial coordinate axis ratio are 
fixed to y4 = 1.0 and r^/re = 0.35, respectively. The circle denotes 
roughly the separation between stable (at the left of the circle) and 
unstable (at the right) models along each sequence. 



FIG. 2: J/M^ of the most rapidly rotating models with a central rest- 
mass density equal to that of the maximum-mass nonrotating models 
{i.e., for the models marked with the circles in Fig.[T) as a function of 
the rotation law parameter A and for different values of A*'. All mod- 
els have J /M^ < 1 indicating the difficulty of constructing unstable 
supra- Kerr models. 



In order to investigate further the effect of the differential- 
rotation-law parameter A and of the EOS on the above conclu- 
sion, we have investigated a large number of rapidly rotating 
models, spanning a wide range of values for A (between 0.6 
and 1.8) and a wide range of polytropic indices (between 0.5 
and 1.5). In all cases, we have computed the value of J/A/^ 
of the most rapidly rotating models we could construct (which 
was normally close to the mass-shedding limit, when it exists) 
for a central rest-mass density equal to that of the maximum- 
mass nonrotating model {i.e. , for the models marked by circles 
in Fig. [Til. The results of this analysis are collected in Fig.|2l 
which shows that all the models with a central density equal to 
the maximum-mass nonrotating stars have J /AP < 1, when 
N and A are allowed to vary. This remains true also when 
considering the unstable models in Fig. [T](;.e., the ones with 
an higher central density), which all have a value of J/M^ 
lower than the models shown in Fig.|2] It is therefore evident 
that no combination of N and A could yield an unstable supra- 
Kerr model. All of these results provide strong evidence that 
all supm-Kerr models examined here are stable. 

Bearing this in mind, it should be noted that because our 
numerical method does not reach exactly the mass-shedding 
limit for any degree of differential rotation (it is difficult to 
achieve convergence at very small values of the axis ratio 
Tp/ Te) and since the existence of a bifurcation between quasi- 



spheroidal and quasitoroidal models"' with the same axis ratio 
and central density has not been investigated yet, we cannot 
strictly exclude the existence of supra-Kerr unstable models. 

The rapidly rotating models with N ~ 1.0, N — 0.75 and 
N = 0.5 shown in Fig.[T]are also shown in Fig.[3]in a mass vs 
maximum energy density plot. Since the most rapidly rotat- 
ing models with differential rotation and small axis ratio are 
quasitoroidal, the corresponding maximum energy density is 
larger than the central energy density by a factor of roughly 
2, depending on the degree of differential rotation. It is not 
yet known whether the value of the central energy density or 
of the off-center maximum energy density is more important 
in determining the stability to axisymmetric perturbations of 
quasitoroidal models. Therefore, the models shown in Figs.|2] 
and [3] could either be only marginally stable or unstable or 
strongly unstable. Nevertheless, the fact that the central den- 
sity of models in Fig. [T]with J/AP > 1 is at least a fac- 
tor of 3 smaller than the central density of the corresponding 
maximum-mass nonrotating models, indicates that even if all 
models in Fig. [3] are well inside the dynamically unstable re- 
gion, there should still be no supra-Kerr unstable models for 
the parameter range examined. 



^ We define as "quasi-spheroidal" those models having the central and max- 
imum rest-mass density being coincident, while we define as "quasi- 
toroidal" those models having the maximum of the rest-mass density not 
located at the center of the stai'. 
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TABLE I: Initial data for the different stellar models. The different columns refer, respectively, to the central rest-mass density pc and its 
maximum pmax, the ratio of the polar to the equatorial coordinate radii Vp/re, the total gravitational mass M, the circumferential equatorial 
radius Re, the central angular velocity Q,c, the ratio of rotational kinetic energy to gravitational binding energy r/| W^|, the ratio J/M^ , where 
J is the angular momentum, and the degree of differential rotation A, where for A — > oo uniform rotation is recovered. All the initial models 
have been computed with a polytropic EOS with K — 100 and N = 1. The last column shows instead the ratio Jbh / M%jj for the BH 
formed after the collapse of the Sub-Kerr models computed using Eq. (5.2) in Ref. (3]. 



Model 


Pc X 10"^ 


/5max X 10"^ 




M/Mq 


Re/M 




T/\W\ 




A 


Jbh/M%u 


Al 


3.0623 


6.8920 


0.23 


1.7626 


3.5424 


0.51891 


0.18989 


0.75004 


0.6 


0.74 


A2 


3.0623 


4.0236 


0.33 


2.2280 


3.5316 


0.21752 


0.21705 


0.81507 


1.0 


0.81 


A3 


3.0623 


3.0623 


0.33 


2.6127 


4.1111 


0.10859 


0.23163 


0.88474 


1.4 


0.88 


Bl 


0.4630 


0.4632 


0.39 


1.9009 


8.8185 


0.03723 


0.21509 


1.08650 


1.0 






FIG. 3: The solid, short-dashed and long-dashed lines (upper curves 
with individual models marked) represent the gravitational mass M 
of some of the unstable models shown in Fig.|2l as a function of the 
maximum energy density, for A'^ = 0.5, 0.75 and 1.0. The different 
values of A are reported near each model. The lower curves show the 
corresponding nonrotating sequence of models for each of the above 
EOSs. 



C. Supra- and sub-Kerr initial data 



FIG. 4: Comparison between our initial models (filled triangles; see 
TableUll and the uniformly rotating models studied in Ref. lil4il . Here 
we plot the gravitational mass M as a function of the central energy 
density Cc. The solid, dashed, and dotted lines correspond to the se- 
quence of nonrotating models, the sequence of models rotating at the 
mass shedding limit, and the sequence of uniformly rotating models 
that are at the onset of the secular instability to axisymmetric pertur- 
bations. Also shown are the secularly (open circles) and dynamically 
unstable (filled circles) initial models used in Ref. |[l4ll (see inset). 



We have investigated the dynamics of differentially rotat- 
ing collapsing compact stars by focusing on three sub-Kerr, 
dynamically unstable models and on one supra-Kerr, dynam- 
ically stable model. For the latter, the collapse was triggered 
through an artificially large pressure depletion of 99%, as was 
done in Ref. UtIi . All models are constructed with the poly- 
tropic EOS with iV = 1 and K = 100. The three dynami- 
cally unstable models are labeled Al to A3 and are shown as 
filled triangles in Fig.|4]while their detailed properties are dis- 
played in Table U The central rest-mass density of the three 
models is chosen to be the same as the central rest-mass den- 



sity of the maximum mass nonrotating model for this EOS. 
The degree of differential rotation varies from A — 0.6 to 
A ~ 1.4. The maximum rest-mass density increases with re- 
spect to the central density, as differential rotation becomes 
stronger {i.e., as the relative length scale A becomes smaller), 
as can be seen by comparing the values in the second and third 
columns of Table H] All three models have comparable values 
of J/AP (0.75 to 0.88), of the ratio of rotational kinetic en- 
ergy to gravitational binding energy T/lVKj (0.19 to 0.23) and 
of the total (gravitational) mass M (1.8 to 2.6), while they 
differ significantly in radius (6.4 to 11) and central angular 
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t=0.000 P =0.000 ms t=0.000 P , =0.000 ms 




X/M X/M 



FIG. 5: Snapshots of the rest-mass density p in the equatorial plane (left panels) and in the xz plane (right panels) for model A2. The contour 
lines are drawn for p = 10~''''^-'^°'^''pmax for j = 0, 1, . . . , 8, where pmax is the maximum of p. After the formation of the AH (dashed line) at 
time tah = 6.71Prot,c, Pmax is the maximum of p at f = tah- Time is normalized to the initial central rotation period of the star, Prot.c = 13M. 
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velocity (0.52 to 0.11). Even though the axisymmetric stabil- 
ity of these models cannot be determined from a turning-point 
method, the numerical simulations performed and discussed 
later on, showed that these models are indeed dynamically un- 
stable and collapse without the need of a pressure depletion. 

Finally, the fourth model studied, model Bl in Table U is 
shown as a empty triangle in Fig. |4] and represents a stable 
supra-Kerr model, with comparable mass and T/|Ty|-ratio as 
models Al to A3, but with much smaller central rest-mass den- 
sity and J /AP = 1.09. The complete set of initial data is 
displayed in Fig. |4] which reports the total mass versus the 
energy density at the center of the star Cc- Also shown as use- 
ful references are the sequence of nonrotating models (solid 
line), the sequence of models rotating at the mass-shedding 
limit (dashed line) and the sequence of models that are at the 
onset of the secular instability to axisymmetric perturbations 
(dotted line). Similarly, indicated respectively with open and 
filled circles, are the secularly (SI — S4) and dynamically un- 
stable (Dl — D4) uniformly rotating models used in the col- 
lapse simulations discussed in Ref. lfl4ll (see Fig. |4] inset). 



Al 

A2 
A3 




0.5 



1.5 



t [ms] 



IV. DYNAMICS OF THE COLLAPSE 

We next discuss the dynamics of the matter during the col- 
lapse of the initial stellar models described in the preceding 
section. All the models were studied with different resolutions 
and using fixed mesh refinement techniques. In the case of the 
sub-Kerr models Al, A2, and A3, up to seven refinement levels 
were used, in order to be able to extract the gravitational-wave 
signal in a region sufficiently distant from the sources. 

The supra-Kerr model Bl, on the other hand, was studied 
using three levels because the dynamics is not limited to a 
small central region of the computational domain (the process 
follows several bounces and subsequent collapses) and so the 
finest refinement level was set to be larger with respect to the 
one used for the sub-Kerr models in order to reduce the nu- 
merical error The outer boundaries were then moved at those 
distances that were computationally affordable, but not suffi- 
ciently far away to allow for gravitational-wave extraction as 
in the sub-Kerr cases. For this model, in fact, we resort to a 
gravitational-wave extraction via the quadrupole formula, as 
discussed in Sec. IV Bl 

An ideal-fluid EOS, P = pe(r - 1), with T = 2 was used 
during the evolution of all the models. 



A. Sub-Kerr Collapse 

All the three sub-Kerr models considered Al, A2, A3 show 
the same qualitative dynamics, with the gravitational collapse 
leading to a central black hole. All of them were evolved in 
equatorial and 7r/2 symmetry {i.e., we considered the region 
{a; > 0, y > 0, z > 0} applying reflection symmetry at z = 
and a rotational symmetry at a; = and y = 0) since they 
did not show the development of nonaxisymmetric instabili- 
ties (when evolved without symmetries at lower resolutions). 



FIG. 6: Maximum of the rest-mass density p normalized to its initial 
value. The triangle, the circle and the square represent the moment 
at which an apparent horizon is formed for models Al, A2, and A3, 
respectively. This figure should be compared with the corresponding 
Fig. [To] which refers to the supra-Kerr model Bl. 



in a way similar to the uniformly rotating models studied in 
Ref. id. 

Because of the similar behavior, we concentrate here only 
on the description of model A2. The results shown here were 
produced on a grid with boundaries located at [0, 137. 9A/] x 
[0, 137. 9M] X [0, 137. 9A/] with a resolution ranging from 
Ax* = 0.86il/ on the coarsest grid to Ax* = 0.013M on 
the finest one, using seven refinement levels. The boundaries 
of the finest grid were chosen in order to include all the star 
and thus to reduce numerical errors. The results shown here 
were obtained without the introduction of any initial perturba- 
tion, except from the truncation error, but a similar dynamics 
was obtained when the collapse was triggered by reducing the 
pressure by 2%, as was done in the case of uniformly rotating 
models in Ref. iflil . 

As one can see in the first row of panels in Fig. |5] where 
we show the isodensity contours in the (x, y) {i.e., equatorial) 
and (a;, z) planes, the star has a toroidal shape due to its strong 
differential rotation. Its evolution is rather similar to what was 
already observed for the uniformly rotating models and espe- 
cially for model D4 in Ref. fl?]. The collapse is axisymmetric 
and leads to the formation of a black hole, as indicated by the 
appearance of an apparent horizon (AH) llssl . The AH, rep- 
resented by a dashed line in Fig.|5] is found at i = 6.71Prot,c, 
where P^oi.c is the initial rotational period at the center of the 
star and is equal to ~ 13i\/. At the time the AH is formed, the 
star has assumed the shape of a disk which rapidly accretes 
until no matter is left outside, as one can see from the last row 
of panels in Fig. |5] Even if we used an ideal-fluid EOS, we 
did not see the formation of global shocks during the collapse. 
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FIG. 7: Total rest-mass A4, and angular momentum J of the matter, 
normalized at t = for the three different models Al (solid line), 
A2 (short-dashed line), and A3 (long-dashed line). The triangle, the 
circle, and the square represent the moment at which an apparent 
horizon is formed for, respectively, model Al, A2, and A3. 



In Fig.|6]we plot instead the maximum of the rest-mass den- 
sity normalized at its initial value for all the three sub-Kerr 
models. All the models show the same dynamics with an ex- 
ponential increase in the maximum of the rest-mass density 
and with the AH forming when p,^^^ has reached a value ~ 2 
times larger than the initial one. As expected, models with a 
higher value of J/AP collapse on a longer timescale. 

In Fig.|2]we compare the total rest-mass and the total angu- 
lar momentum of all the three models (Al, A2, A3) normalized 
to their initial values. As one can see from this figure, the level 
of accuracy for these simulations is very high, with a relative 
error in the conservation of rest-mass and angular momentum 
smaller than 10^'^. Note also that the angular momentum is 
computed from the matter sources, outside the apparent hori- 
zon, and hence it is shown to vanish at late times in Fig. |71 
as the matter is dissipated near the singularity (see discussion 

in illMi). 



B. Supra-Kerr Collapse 

Model Bl has J/M^ ~ 1.1 and shows very different dy- 
namics with respect to the sub-Kerr models. Because in this 
case we expected the development of nonaxisymmetric in- 
stabilities, we decided to adopt equatorial and 7r-symmetry 
(this means that we evolved only the region {x > Q,z > 
0} applying a rotational symmetry boundary condition at 
x = and reflection symmetry at z = 0). Here we re- 
port the results obtained on a grid with boundaries located at 
[0,34A/] X [-34A/,34A/] x [0,34A/] and with a resolution 



ranging from Ax' = 0.17M to Aa;' = OMM with three 
refinement levels and with the finest grid covering the entire 
star. It is worth remarking that the resolutions adopted here for 
the supra-Kerr model are considerably finer than the ones rou- 
tinely adopted in the simulation of binary neutron stars (see, 
e.g., iSliSinillTtl), andof a factor ~ 4 larger As a result, 
although this scenario is much simpler to simulate, its compu- 
tational costs are indeed equally high. 

Because model Bl is a very stable configuration, we had 
to enforce its collapse by artificially reducing the initial pres- 
sure by 99%, as was done by Duez et al. lUTT . while smaller 
pressure reductions were found to be insufficient to trigger 
the collapse. With its pressure support removed, the model 
immediately flattens along the z-direction and collapses to- 
ward the center on the equatorial plane, producing a strong 
shock. After a first bounce, due to the centrifugal barrier pro- 
duced by the large angular momentum, a quasitoroidal struc- 
ture forms, which rapidly fragments into four clumps (see the 
snapshots from time t = 1.122Pi-ot,c to t = 2.154R 

ot,c m 

Fig. [8]) whose formation was observed also in Ref. lUTIl . We 
have also extracted the Fourier modes of the rest-mass den- 
sity p, by computing at z = and at different cylindrical radii 

137 = x'^ + j/^, the averages 



p{w' cos(0), w' sin(0))e""'^d0 . (13) 



The mode power P,„ is then simply given by 
1 



^out — ^in Jro: 



\k^{w)\dw, (14) 



where zuin and t^out are chosen to cover the whole domain 
(for details and an extensive use of this technique, see also 
Refs. iMUl). 

The presence of an m = 4 mode at the beginning can be 
then seen looking also at the modes' power plotted in Fig. |9] 
where we show the evolution of the m ~ 2 (solid line) and 
?7i = 4 (dashed line) modes, being the other modes either 
zero or much smaller during the entire simulation. It is not 
clear at the moment whether the fragmentation has to be con- 
sidered physical and only triggered by the use of a Cartesian 
grid, or entirely due to our Cartesian coordinate system. We 
note that Truelove and collaborators lIHoll have shown that spu- 
rious fragmentation can occur if the Jeans length is not well 
resolved, i.e. if the following "Jeans condition" is verified. 



Ax ^ 1 
a7 ^ 4' 



where Xj is the Jeans length and is given by 




(15) 



(16) 



and Cs is the sound speed. Duez et al. lllTll estimated the mini- 
mum of the Jeans length to be Aj « 1.3M for a model similar 
to our model Bl and using a poly tropic EOS (for an ideal- 
fluid EOS, as the one used in our simulations, the sound speed 
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FIG. 8: Snapshots of the rest-mass density p in the equatorial xy plane (left panels) and yz plane (right panels) for model Bl. The contour 
lines are drawn for p = 10~'''''^-'^'''^'max(p) for j — 0, 1, ... , 12. Time is normalized to the initial central rotational period of the star, 
Protx = SMI. 
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FIG. 9: Logarithm of the power in the m — 2 (solid line) and m = 4 
(dashed line) Fourier modes as a function of time for model Bl. The 
modes not shown in this plot are zero during all of the evolution. 

is generically larger). In the simulation of Duez et al. iIitI . 
but also in ours, the value of Ax/Xj is indeed found to be 
smaller than 0.25, leading to the conclusion that the fragmen- 
tation is physical and due to a genuine nonaxisymmetric in- 
stability. We should remark, however, that in Ref. ||60ll the 
condition ( fTSl l was indicated as necessary but not in general 
sufficient to avoid the formation of spurious fragmentation. 
Thus, even if the resolutions used in our simulation and in the 
one reported in Ref. 1 17] satisfy the Jeans condition, we can- 
not thus strictly exclude that the origin of this m = 4 mode 
is due to the use of a Cartesian grid. Further investigations 
with resolutions higher than the ones that could be afforded 
here and with different coordinate systems will be necessary 
in order to clarify this issue. 

At time t 2.5Prot,c the four fragments merge and a new 
collapse and bounce follows with the formation of a new qua- 
sitoroidal structure. The effects of these bounces on the max- 
imum of the rest-mass density p and specific internal energy e 
are shown in Fig. [TO] At i « 3.0Piot,c the quasitoroidal struc- 
ture contracts toward the center, forming a new configuration 
which does not collapse further, but that develops a bar which 
lasts for w 2 ms. At this point, the model approaches a new 
stable configuration, as one can easily see from the worldline 
of the maximum of the rest-mass density p (Fig.fTTTi and from 
its evolution (Fig.fTOl top panel). It is also evident from Fig.|9] 
that at late times the m = 2 mode dominates. 

In Fig.[T2]we plot the evolution of the total rest-mass and of 
the total angular momentum normalized to their initial values. 
When we stopped the simulation, J had dropped by ~ 10%. 
The loss in the angular momentum cannot be accounted for in 
full by the emission of gravitational waves and it is due to the 
loss of mass, which is expelled by shocks through the external 



FIG. 10: The maximum value of rest-mass density p and of the spe- 
cific internal e normalized to their initial value and during the evo- 
lution of the supra- Ken' model Bl. Note that the end-configuration 
is that of equilibrium, with larger central density and internal energy. 
This figure should be compared with the corresponding Fig.|6] which 
refers to the sub-Kerr models Al — A3. 

boundaries. In the inset of Fig. [12] we plot instead the profiles 
of the angular velocity along the x-axis, fl ~ {av^ — 
at the initial time (purple long-dashed line) and at the end 
of the simulation (light blue solid line), and where is the 
three-velocity of the fluid as measured by an Eulerian ob- 
server lITill . For comparison we also show the Keplerian pro- 
file oc r~'^/^ (black short-dashed line). It is evident that 
the collapse leads to an even larger degree of differential ro- 
tation than the initial one and that the new equilibrium is very 
close to a Keplerian configuration in its outer layers. This 
is clearly the result of having essentially removed pressure 
forces, leaving the centrifugal ones the only ones responsible 
for the equilibrium. 

As a final remark we note that the fact that we were not 
able to force this model to collapse to a black hole, even 
when artificially reducing the pressure by 99%, confirms that 
supra-Kerr models cannot directly collapse to a black hole. 
Furthermore, the evidence that even when forced to collapse 
the supra-Kerr model does not produce a black hole, but 
rather redistributes its angular momentum to reach a stable 
and axisymmetric stellar configuration, provides strong evi- 
dence that cosmic censorship is not violated and that rather 
genetic conditions for a supra-Kerr progenitor do not lead to 
a naked singularity. 

V. GRAVITATIONAL- WAVE EMISSION 

We now concentrate on the emission of gravitational waves 
from the sub-Kerr and supra-Kerr models, with the aim of 
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FIG. 1 1 : Worldline of the maximum of the rest-mass density p for 
model Bl. The final configuration has a maximum located at a radius 
smaller than the initial one. Note also the three bounces. 

comparing our results with those obtained in Refs. iIztI l6lll 
for the collapse of uniformly rotating neutron stars models. 

A. Sub-Kerr Models 

In Fig.[T3]we show hj^ and for the A2 sub-Kerr model. 
We recall that in the case of axisymmetric collapse to black 
hole (BH) the signal is expected to be composed by an ini- 
tial rapid increase in the amplitude, due to the initial phase 
of the co llap se, followed by a ring-down phase of the formed 
BH IztI [35I I6ll - l64ll . These two phases are very evident in 
Fig- Eland the same behavior is shown also by the other sub- 
Kerr models. The initial oscillations visible in the /i+ for 
the first 0.5ms are spurious effects related to the initial vio- 
lation of the Hamiltonian constraint introduced by importing 
axisymmetric models computed on a spherical grid onto the 
Cartesian grid used for the evolution. We have verified indeed 
that these oscillations decrease in amplitude with increasing 
resolution, while the rest of the signal remains the same. 

In Table we report the signal-to-noise ratio for various 
second and third-generation gravitational-wave detectors and 
the energy emitted through gravitational waves for all the 
sub-Kerr models, assuming a galactic source at a distance of 
10 kpc. For example, gravitational waves from the collapse 
of model A2, with a total energy of about 1.4 X IQ-Hl, 
would arrive at the Virgo detector with a signal-to-noise 
ratio of S/N w 2.11 (with a characteristic amplitude of 
he = 6.14 X 10^^^(A/ /Mq) and at a characteristic frequency 
fc = 1367Hz). In the case of LIGO, instead, we obtain he = 
5.28xlO~2i(M/AfQ)at/c = 1133Hz with 5/iV « 1.5 while 
for advanced LIGO we have he = 3.81 x IQ-^^M/Mq) 



FIG. 12: Total rest-mass M, and total angular momentum J nor- 
malized at their initial values for the supra-Kerr model Bl. The in- 
set shows instead the profiles of the angular velocity at the initial 
time (purple long-dashed line) and at the end of the simulation (light 
blue solid line). For comparison we also show the Keplerian profile 
oc r""^/^ (black short-dashed line). 

at fc ^ 809Hz with S/N w 14.7. However, when con- 
sider third-generation detectors such as the Einstein Telescope 
(ET) ll65ll . the signal-to-noise ratio would increase dramati- 
cally to S/N 247 and the source could then be detected up 
to a distance of ~ 1 Mpc. Even though detection of individual 
events is difficult, the case of a stochastic background is worth 
further investigation. 

The gravitational-wave energies computed here are simi- 
lar to those obtained for the collapse of uniformly rotating 
neutron stars, as reported in ll6l1l . This is to be expected, 
since even though the initial quadrupole moment is larger 
for differentially rotating models, the collapse proceeds more 
slowly, due to the higher values of J/AP and this reduces 
the efficiency in energy emission through gravitational waves. 
This can also be seen in Fig. [141 where we plot the total en- 
ergy emitted by gravitational waves, AE/M, as a function 
of J'^/M, which is proportional to the initial axisymmetric 
quadrupole moment. The triangles represent the uniformly ro- 
tating neutron stars discussed in |6l]], the squares represent the 
differentially rotating models discussed here, while the solid 
line is the best fit of the data with the following analytic ex- 
pression: 

AE (JVM)"- 

where ai = (5.17 ± 4.37) x 10^ 03 = (1.11 ± 0.57) x 10^ 
?i2 = 2.63 ± 0.53, and ni = 1.43 ± 0.74. Clearly, because 
of the small statistics, the error in some of the coefficients 
is rather large, but this can be compensated by performing 
additional simulations. 
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TABLE II: Signal-to-noise ratio (SNR) computed for the collapse of sub-Kerr and supra-Kerr models assuming a source at a distance of 
10 kpc. The last column shows also the energy emitted through gravitational waves which, in the case of the sub-Kerr models, is comparable 
to the one obtained in (KUl from the collapse of uniformly rotating neutron stars. 



Model 


SNR (Virgo) 


SNR (LIGO) 


SNR (Adv. Virgo) 


SNR (Adv. LIGO) 


SNR (ET) 


AE/M 


Al 


1.00 


0.60 


3.56 


4.80 


107.41 


5.46 X lO"'^ 


A2 


2.11 


1.47 


13.98 


14.67 


246.77 


1.38 X lO"'' 
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4.55 


3.79 


45.70 
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594.54 


6.24 X 10~* 
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16.13 
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FIG. 13: Gravitational-wave amplitudes h+ (top panel) and hx (bot- 
tom panel) for the collapse of model A2. The waves are computed 
using gauge-invariant perturbations of a Schwarzschild spacetime. 



B. Supra-Kerr Model 



FIG. 14: Energy emitted through gravitational waves, normalized to 
the total initial mass AI, as a function of /M. Triangles represent 
the uniformly rotating neutron star models studied in ll6lll . while the 
squares refer to the differentially rotating models discussed here. The 
solid line is the best fit (see text for details). 



As anticipated in Sec. |IV] we have also estimated the 
gravitational- wave signal emitted by model Bl using the stan- 
dard quadrupole formula ll66l l67ll since in this case the outer 
boundary was not located far enough from the source to allow 
for the extraction of the signal using the method described in 
Sec. 



IIIDI In this approximation, the observed waveform and 
amplitude for the two polarizations for an observer situated 
at large distance r along the z-axis are approximately given 
by Q 



hx = 2 



r 

ry{t) 



where 



= (PxD x^x" 



(18) 
(19) 

(20) 



is the quadrupole moment of the matter distribution, D = 
y/jpW, 7 is the determinant of the three metric, and W is 
the Lorentz factor The results for the gravitational-wave am- 
plitudes /i+ and hx are finally reported in Fig. [15] and they 
are 1 order of magnitude larger than those emitted by the sub- 
Kerr models. This produces consequently much larger signal- 
to-noise ratios as one can see by looking at the last row of 
table ini As an example, if such a collapse happened at a 
distance of 10 kpc, we would obtain in the case of advanced 
LIGO a signal-to-noise ratio of S/N « 72, with a charac- 
teristic amplitude of he = 1.32 x 10"^^(M/Afo) and at a 
characteristic frequency fc = 2714Hz. In the case of ET in- 
stead the signal could be detected for a source at a distance of 
up to w 10 Mpc. 
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FIG. 15: Gravitational- wave amplitudes h+ and hx for the col- 
lapse of model Bl. The waves are computed using the Newtonian 
quadrupole formula. 

VI. CONCLUSIONS 

We have presented new results on the collapse of differ- 
entially rotating neutron stars in full general relativity, using 
simulations with fixed mesh refinement and in 3D. We have 
considered two different classes of initial configurations, con- 
sisting of either sub-Kerr models with J/A/^ < l,orofsupra- 
Kerr models with J /AP > 1. In particular, we have per- 
formed a detailed study of the parameter space of different 
equilibrium configurations and, for polytropes which could 
represent neutron-star models, we could only find models un- 
stable to an axisymmetric instability that are sub-Kerr. 

When evolving three representative sub-Kerr models with 
different degrees of differential rotation we observed in all 
cases a dynamics which is very similar to the one already 
found in the case of uniform rotation. In particular, all mod- 
els remain axisymmetric during collapse and produce a Kerr 
black hole. We have also studied the collapse of a supra-Kerr 
model when its collapse is triggered by a 99% depletion of its 
pressure. In this case, we observed a very different dynamics, 
with the formation of nonaxisymmetric instabilities, the de- 
velopment of a torus and its subsequent fragmentation in four 
clumps that merge again, forming a bar and eventually a stable 
axisymmetric configuration. While lack of evidence cannot be 



taken to exclude the possibility that a naked singularity can be 
produced by the collapse of a differentially rotating star, our 
results also suggest that cosmic censorship is not violated and 
that rather generic conditions for a supra-Kerr progenitor do 
not lead to a naked singularity. 

For all the models presented here, we have also computed 
the gravitational-wave signal, which is comparable to the one 
from uniformly rotating neutron stars. In particular, the effi- 
ciency in the emission of gravitational waves for the sub-Kerr 
models is slightly smaller than the one of uniformly rotating 
stars. This is due to the large angular momentum that differ- 
entially rotating neutron stars have and which increases the 
time scale on which the collapse happens. Therefore, even 
if the initial quadrupole moment is much larger, the longer 
time needed for the collapse decreases the overall efficiency. 
Furthermore, as for the collapse of uniformly rotating neutron 
stars, also in this case the detection of such events would be 
possible only for sources located in our Galaxy in the case of 
the advanced Virgo and advanced LIGO detectors. However, 
a third-generation detector such as ET could instead be able 
to detect a source at a distance of ~ IMpc. 

We plan to extend the work presented here by making use 
of the axisymmetric code developed in Ref. ll68ll . which will 
allow us to consider a larger number of models and EOSs. We 
also plan to investigate the effect that magnetic fields have on 
these configurations and especially on the supra-Kerr models 
by using the GRMHD version of Whisky 12811 . The magnetic 
field can in fact redistribute the angular momentum inside the 
star and thus alter the dynamics of the collapse, as already 
studied in axisymmetry ll20l422ll and in 3D llioll . 
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